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Abstract.  Carbon  related  defects  are  readily  incorporated  in  GaN  due  to  its  abundance 
during  growth  both  with  MBE  and  CVD  techniques.  Employing  first-principles  calculations 
we  compute  the  migration  barrier  of  neutral  carbon  interstitials  in  the  wurtzite  GaN  crystal. 
The  Minimum  Energy  Path  (MEP)  and  the  migration  barriers  of  these  defects  are  obtained 
using  the  Nudged  Elastic  Band  (NEB)  method  with  the  climbing  image  modification  (CI-NEB). 
In  addition,  the  Dimer  method  is  used  to  verify  the  results.  The  results  yield  a  quantitative 
description  of  carbon  diffusion  in  the  crystal  allowing  for  the  determination  of  the  most  probable 
migration  paths. 


1  Introduction 

Efficient  LED  lighting  has  developed  in  recent  years  thanks  to  Ill/nitride  semiconductors  such 
as  GaN  [1].  In  addition  to  optoelectronic  applications  such  as  blue/green  LEDs,  GaN  and  its 
related  alloys  are  used  in  power  electronics  [2]  as  well  as  photovoltaic  applications  [3]. 

Defects,  both  point  [4]  and  extended  ones  [5,  6],  affect  the  electronic  properties  of  devices. 
Among  them,  carbon  is  a  relevant  and  easily  incorporated  defect  in  GaN.  Defect  related  levels 
in  the  band-gap  may  be  the  source  of  radiative  recombination  centers,  leading  to  below  gap 
emission  compromising  the  performance  of  the  device.  A  typical  example  of  such  a  level  is 
centered  around  2.2  —  2.3  eV  and  is  often  referred  to  as  the  yellow  luminescence  (YL)  [7-10]. 
This  YL  band  can  be  present  for  both  undoped  samples  [11]  and  samples  containing  carbon 
impurities  [12-14].  Recent  studies  suggest  that  carbon  related  defects  are  responsible  for  the 
YL  band  [7,  15]. 

Carbon  is  a  common  impurity  both  in  molecular  beam  epitaxy  (MBE)  [16,  17]  and  metal 
organic  chemical  vapor  deposition  (MOCVD)  [18].  In  the  former  case,  carbon  can  contaminate 
the  material  during  air  exposure  in  standard  substrate  loading  procedures  or  at  the  beginning  of 
regrowth.  In  MOCVD,  carbon  is  part  of  the  metal  organic  compounds  used  as  source  material 
for  gallium.  In  addition,  carbon  can  be  found  as  a  contaminant  in  the  source  gases  or  it  can  be 
etched  off  the  susceptor  that  transfers  heat  to  the  substrate.  Even  though  carbon  defects  have 
been  intensively  investigated  in  the  past,  the  current  understanding  of  the  migration  mechanisms 
of  these  defects  in  GaN  is  still  incomplete. 

The  migration  paths  and  barriers  of  carbon  interstitials  in  GaN  were  investigated  using 
density-functional  theory  (DFT)  [19,  20]  in  conjunction  with  the  climbing  image  nudged  elastic 
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band  method  (CI-NEB)  [21,  22].  Using  the  CI-NEB  method  we  ensure  that  the  calculated 
energy  corresponds  to  the  migration  barrier  for  a  given  path.  In  addition,  we  get  information 
about  the  atomic  configuration  in  the  saddle  point. 

2  Background 

It  is  empirically  observed  that  most  diffusion  processes  in  solids  exhibit  a  temperature 
dependence  described  by  an  Arrhenius  law 


D  =  Doe  k*T,  (1) 

where  Dq  is  a  temperature  independent  factor,  Q  is  the  activation  energy  for  the  atomic  jump 
mechanism  and  k^T  is  the  product  of  the  Boltzmann  constant  with  temperature. 

Deviations  from  this  equation  may  be  caused  by  the  existence  of  short  circuiting  diffusion 
paths  such  as  dislocations  and  grain  boundaries,  multiple  diffusion  mechanisms  and  impurity 
effects  [23]. 

The  jump  rate  is  characterized  by  an  attempt  frequency  and  a  thermodynamic  factor  that 
dictates  the  probability  of  an  attempt  resulting  in  a  successful  jump. 

AG 

r  =  r0e  k*T.  (2) 

The  attempt  frequency,  Tq,  can  be  obtained  within  harmonic  transition  state  theory  via  the 
Vineyard  equation  [24],  but  often  the  Einstein  or  Debye  frequency  is  used  instead.  AG  represents 
the  isothermal  work  associated  with  the  motion  of  the  particle  from  the  initial  to  the  final  state. 

3  Computational  Details 

3.1  Computational  method 

The  calculations  were  performed  with  the  Vienna  ab  initio  simulation  package  (VASP)  [25]  using 
the  projector  augmented  wave  (PAW)  method  [26,  27].  The  generalized  gradient  approximation 
(GGA)  in  the  parametrization  by  Perdew,  Burke  and  Ernzerhof  (PBE)  [28]  was  used.  The 
investigation  of  the  lowest  energy  paths  was  performed  using  the  CI-NEB  method  [21,  22]  as 
well  as  the  Dimer  method  [29]  as  implemented  in  VASP  through  the  VTST-Tools  by  Henkelman, 
Jonsson  and  others  [30]. 

The  convergence  of  the  size  of  the  supercell  was  investigated  using  supercells  of  64,  72,  96,  128 
and  144  atoms.  The  64-atom  supercell  is  found  to  be  inadequate  to  contain  the  deformations 
induced  by  the  migration  of  carbon.  The  72-atom  supercell  might  be  sufficient  for  certain 
migrations  but  not  all.  Thus,  we  chose  the  96-atom  supercell  to  be  the  safest  and  most  efficient 
option  at  hand.  The  plane-wave  energy  cutoff  was  set  at  450  eV  and  a  T-centered  2x2x2 
k-point  mesh  is  used  for  the  sampling  of  the  Brillouin  zone.  The  atomic  configurations  were 
relaxed  until  the  maximum  force  per  atom  was  less  than  5  x  10-3  eV /A.  In  the  case  of  the 
NEB  calculations,  the  images  were  relaxed  until  the  maximum  force  per  atom  was  no  more  than 
10-2  eV/A. 

The  crystallographic  parameters  of  the  w-GaN  we  obtained  from  our  calculations  are 
a  =  3.211  A,  c/a  =  1.629  and  u  =  0.377,  which  are  in  good  agreement  with  the  experimental 
values  a  —  3.189  A,  c/a  —  1.626  and  u  —  0.377  and  consistent  with  previous  DFT  calculations 
[31].  The  bandgap  at  the  T  point  is  1.757  eV  which  underestimates  the  experimental  value  of 
3.4  eV  as  expected  from  GGA. 

The  main  sources  of  error  in  DFT  calculations  of  point  defects  are  the  electrostatic  and 
elastic  interactions  of  the  defect  between  neighboring  supercells  and  the  underestimation  of  the 
bandgap.  Lany  and  Zunger  [32]  describe  a  number  of  methods  to  overcome  these  shortcomings 
of  DFT.  However,  migration  barriers,  the  main  focus  of  the  present  work,  are  calculated  as 
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(a)  C-N  type-1  split  interstitial  (b)  C-N  type-2  split  interstitial 


Figure  1:  Split  interstitials  of  carbon  (yellow)  and  nitrogen  (blue)  surrounded  by  four  gallium 
atoms  (red). 


energy  differences  of  electronically  similar  configurations.  Thus,  there  is  no  need  to  apply  any 
correction  scheme. 

A  typical  systematic  error  in  the  formation  energies  calculations  is  considered  to  be 
approximately  0.1  eV  [33].  Performing  identical  migrations  at  different  sites  we  verified  that 
the  uncertainty  of  the  migration  barrier  also  lies  within  0.1  eV. 

In  the  NEB  method  [34],  a  set  of  “images”  of  the  system  is  used  to  represent  the  migration 
path  from  the  initial  to  the  final  configuration.  The  images  (i.e.,  atomic  configurations  along 
the  migration  path)  are  connected  with  springs  to  resemble  a  string  (or  band).  An  optimization 
algorithm  is  then  applied  to  relax  the  string  down  towards  the  minimum  energy  path.  In 
the  climbing  image  modification  the  highest  energy  image  is  driven  to  the  saddle  point  by 
neutralizing  the  forces  along  the  band.  Since  the  image  converges  to  the  saddle  point,  the  exact 
migration  barrier  can  be  calculated.  In  addition,  the  Dimer  method  was  used  to  verify  the  saddle 
point  location  and  the  migration  barrier  acquired  via  the  CI-NEB. 

3.2  Determination  of  diffusion  paths 

The  two  most  stable  configurations  for  the  neutral  carbon  interstitial  are  the  type-1  and  type-2 
C-N  split  complexes  shown  in  Fig.  1.  They  differ  by  less  than  0.1  eV  [35].  Hence,  both  cases 
are  investigated  as  potential  initial  and  final  configurations  for  the  diffusion  of  carbon. 

This  small  energy  difference  and  the  similarity  of  the  two  configurations  allows  for  migrations 
accompanied  by  transformation  from  one  type  to  the  other.  In  this  case,  the  migration  barrier 
is  defined  as  the  energy  difference  between  the  saddle  point  and  the  lowest  energy  configuration. 
As  a  result,  the  barrier  in  this  case  could  be  smaller  by  ^  0.1  eV  following  the  opposite  direction 
of  the  reaction. 

The  migration  paths  investigated  in  the  present  work  have  components  both  parallel  and 
perpendicular  to  the  c-axis  of  the  wurtzite  crystal.  The  endpoints  of  the  migrations  are  either 
a  type-1  or  type-2  split  interstitial.  Fig.  2  presents  the  different  migration  paths  considering 
jumps  among  the  first  and  second  nearest  neighbors.  The  different  combinations  of  paths  and 
endpoint  configurations  result  in  nine  different  calculations. 

Process  A,  shown  in  Fig.  2  refers  to  motion  between  two  out  of  plane  first  nearest  neighbors. 
A  single  gallium  atom  forms  four  bonds  with  nitrogens  in  the  ideal  crystal.  Three  of  these 
nitrogens  lie  on  a  plane  which  is  perpendicular  to  the  [0001]  direction.  Process  A  describes 
the  migration  of  a  carbon  atom  forming  a  split  with  one  of  these  three  nitrogens  to  the  fourth 
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C  -  2nd  nbr  out  of  plane 
A  -  1st  nbr  out  of  plane 


B  -  1st  nbr  in  plane 


Figure  2:  Diffusion  paths  available  to  carbon  interstitial  considering  jumps  among  the  first  and 
second  nearest  neighbors. 


nitrogen  which  is  out  of  plane. 

In  the  case  of  process  B,  migration  occurs  between  two  adjacent  nitrogens  lying  on  a  plane 
perpendicular  to  the  [0001]  direction.  This  process  is  a  first  nearest  neighbor  migration  process 
and  together  with  process  A  is  expected  to  exhibit  the  lowest  migration  barriers. 

As  an  additional  process,  we  consider  process  C  in  which  migration  occurs  between  the  second 
nearest  neighbors  which  lie  in  different  planes  with  respect  to  the  [0001]  direction.  Even  though 
this  process  utilizes  the  hexagonal  channel  of  the  wurtzite  structure,  it  is  expected  to  exhibit 
the  highest  barriers  since  it  is  a  second  nearest  neighbor  migration. 

4  Results 

The  difference  of  the  formation  energy  between  the  two  types  of  interstitials  is  0.13  eV  with  the 
type-2  being  the  lowest  in  energy.  The  bond  length  between  the  carbon  and  the  nitrogen  is 
1.29  A  and  1.31  A  for  the  type-1  and  type-2  respectively. 

The  results  obtained  using  the  CI-NEB  method  are  summarized  in  Table  1.  The  lowest 
migration  barrier  is  observed  in  the  case  of  process  A  from  a  type-2  interstitial  to  a  type-1. 
This  barrier  becomes  even  smaller  by  approximately  0.1  eV  if  the  reaction  path  is  reversed.  The 
diffusion  potential  energy  along  the  reaction  path  is  shown  in  Fig.  3. 

The  Dimer  method  verifies  the  results  acquired  by  the  CI-NEB  for  all  the  process  A 
migrations.  Both  methods  result  in  the  same  atomic  configuration  for  the  saddle  point  and 
consequently  the  same  migration  barrier.  In  the  case  of  the  lowest  barrier  migration,  the  saddle 
point  lies  closer  to  the  type-1  interstitial.  With  respect  to  the  straight  line  between  the  initial 
and  final  configurations,  the  saddle  point  forms  an  angle  of  29°  and  16°  with  the  type-1  and 
type-2  ends  respectively. 

First  nearest  neighbor  in  plane  migrations  (process  B)  exhibit  barriers  close  to  3eV  making 
them  unfavorable  compared  to  their  out  of  plane  counterparts.  As  expected,  second  nearest 
neighbor  out  of  plane  migrations  require  a  much  greater  activation  energy. 
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Figure  3:  Diffusion  energy  barrier  for  a  process  A  migration  from  a  type-1  to  type-2  split  carbon 
interstitial. 

Table  1:  Energy  barriers  for  carbon  interstitial  migration  as  obtained  from  CI-NEB  calculations. 


Migration  Process 

Initial 

Final 

Barrier  (eV) 

si 

si 

2.6 

A 

s2 

s2 

2.7 

si 

s2 

2.3 

si 

si 

2.9 

B 

s2 

s2 

3.1 

si 

s2 

3.1 

si 

si 

4.0 

C 

s2 

s2 

— 

si 

s2 

4.0 

In  order  to  estimate  the  temperature  where  carbon  interstitials  are  activated  we  assume  a 
jump  rate  of  1  Hz  and  a  typical  Debye  frequency  of  lOTHz.  Then,  using  eq.  (2)  and  the 
calculated  migration  barriers,  one  can  estimate  the  temperature  where  migration  occurs.  Hence, 
using  T  =  1  Hz,  Tq  =  1013  Hz  and  AG  —  2.3  eV  we  estimate  that  carbon  interstitials  are  activated 
at  temperatures  close  to  890  K.  This  temperature  is  relevant  during  growth  of  GaN  both  with 
MBE  and  CVD  techniques. 

5  Conclusions 

We  employed  the  Climbing  Image  Nudged  Elastic  Band  technique  to  explore  migration  barriers 
of  carbon  diffusion  in  GaN.  The  Dimer  method  was  also  used  in  order  to  verify  the  results 
obtained  by  the  CI-NEB  in  the  case  of  the  three  lowest  energy  barrier  migrations.  The  two  most 
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stable  interstitial  types  of  carbon  are  investigated  with  respect  to  diffusion  via  first  and  second 
nearest  neighbor  mechanisms. 

The  first  out  of  plane  nearest  neighbor  mechanism  exhibits  the  lowest  migration  barriers 
indicating  a  favorable  direction  of  migration  along  the  [0001]  axis  of  the  crystal.  However, 
for  growth  temperatures,  where  carbon  migration  is  relevant,  the  in  plane  migration  is  also 
competing.  Migration  via  the  second  nearest  neighbor  mechanism  exhibits  high  energy  barriers 
minimizing  the  probability  for  this  mechanism  to  contribute  to  the  carbon  diffusion. 
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